class: inverse, middle # Intro to R - 2 .font100[ Bon Woo Koo & Subhro Guhathakurta 9/1/2022 ] --- ## Content * Rpubs & R Markdown * Loops and functionals * geospatial data in R --- class: inverse, center, middle # Rpubs & R Markdown --- ## Rpubs RPubs is an open platform by R Studio (the company) for publishing R Markdown documents and presentations. * It's free. * Everything you upload there is open to public. * **DO NOT** upload anything you need to protect (e.g., API keys, personal info etc.). To create an account, [**click here**](https://rpubs.com/). --- ## R Markdown **Markdown** is a markup language that adds formatting elements to plaintext text documents. **R Markdown** is a markdown that contains chunks of embedded R code. * Our syllabus is a R Markdown document in which all tables are written in R code. <br> * A cooler example is [this](http://timelyportfolio.github.io/rCharts_nyt_home_price/). To create a R Markdown file, go to **file > R Markdown.**  --- --- ## Publishing on Rpubs  For your assignment, you can submit the link, which in this case would be: http://rpubs.com/BonWooKoo/this_will_your_link --- class: inverse, middle, center # Loops --- ## Loops Loops are used to repeat a specific task. Ones that I commonly use include for loop, while loop, and apply. .pull-left[ **A for loop** iterates over a vector. ```r for (i in 1:3){ print(i) } ``` ``` ## [1] 1 ## [1] 2 ## [1] 3 ``` ```r for (i in 1:3){ print(i + 2) } ``` ``` ## [1] 3 ## [1] 4 ## [1] 5 ``` ] .pull-right[ **A while loop** keeps repeating the task as long as a specified condition is true. ```r n <- 0 while (n < 5){ # if (n >= 0), it will run forever print(n) n <- n + 1 } ``` ``` ## [1] 0 ## [1] 1 ## [1] 2 ## [1] 3 ## [1] 4 ``` ] --- ## Functional Hadley Wickham strongly advocates for *NOT* using loops and advocates functionals, such as **lapply()** or **map()**. .pull-left[ ```r # Loop for (i in 1:3){ print(i+1) } ``` ``` ## [1] 2 ## [1] 3 ## [1] 4 ``` ] .pull-right[ ```r # Functional add_one <- function(x) x+1 map_dbl(1:3, add_one) ``` ``` ## [1] 2 3 4 ``` ] In this course, however, I will try to use loops than functionals. I personally think loops are more intuitive for beginners. If you want to learn more about it, you can read [this](https://adv-r.hadley.nz/fp.html). --- class: inverse, middle, center # Geospatial data in R --- class: center, middle  --- class: center, middle  --- class: center, middle  --- .font180[.green[sf package in R]] The **sf** stands for simple features. **"Simple features are a standardized way of encoding spatial vector data (points, lines, polygons) in computers" (Pebesma, 2018, 439).** .font80[ * *Feature* can be thought of as "things" or objects that have a spatial location or extent (e.g., building or political state). * *Feature geometry* refer to the spatial properties (location or extent) of a feature. * *Feature attributes* refer to other properties that features have, such as name, some measured quantity, etc. ] .footnote[ This content is adopted from Pedesma (2018) [Read](https://journal.r-project.org/archive/2018/RJ-2018-009/RJ-2018-009.pdf) ] --- class: center, middle  --- class: center, middle  --- class: center, middle  --- class: center, middle  --- * Fast reading and writing of data (but not necessarily plotting) * sf objects can be treated as data.frames * Works smoothly with dplyr verbs as well as %>%. * I personally do not use ArcGIS at all. 99% of my GIS works are done in R using sf package.  .footnote[image source: https://r-spatial.github.io/sf/articles/sf1.html] --- ## Let's do some coding ```r library(sf) census <- sf::st_read("https://raw.githubusercontent.com/BonwooKoo/UrbanAnalytics2022/main/Lab/module_0/testdata.geojson") marta <- sf::st_read("https://raw.githubusercontent.com/BonwooKoo/UrbanAnalytics2022/main/Lab/module_0/MARTA_Routes.geojson") five_point <- sf::st_read("https://raw.githubusercontent.com/BonwooKoo/UrbanAnalytics2022/main/Lab/module_0/five_point.geojson") ``` .footnotesize[.scroll-box-18[ ```r census %>% print() ``` ``` ## Simple feature collection with 519 features and 8 fields ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -84.85071 ymin: 33.35246 xmax: -84.02371 ymax: 34.18629 ## Geodetic CRS: WGS 84 ## First 10 features: ## GEOID tract county state hhincome race.tot race.white race.black geometry ## 1 13063040307 Census Tract 403.07 Clayton County Georgia 37064 4939 1442 1784 MULTIPOLYGON (((-84.37445 3... ## 2 13063040512 Census Tract 405.12 Clayton County Georgia 36414 5478 351 4996 MULTIPOLYGON (((-84.44869 3... ## 3 13063040519 Census Tract 405.19 Clayton County Georgia 32568 4380 35 4178 MULTIPOLYGON (((-84.42895 3... ## 4 13063040617 Census Tract 406.17 Clayton County Georgia 33182 1418 941 147 MULTIPOLYGON (((-84.27173 3... ## 5 13067030224 Census Tract 302.24 Cobb County Georgia 73061 6318 4019 1262 MULTIPOLYGON (((-84.65095 3... ## 6 13067030235 Census Tract 302.35 Cobb County Georgia 104265 4862 3904 682 MULTIPOLYGON (((-84.70304 3... ## 7 13067030313 Census Tract 303.13 Cobb County Georgia 93817 9073 6391 1416 MULTIPOLYGON (((-84.55734 3... ## 8 13067030412 Census Tract 304.12 Cobb County Georgia 38668 5018 2221 2165 MULTIPOLYGON (((-84.5101 33... ## 9 13067031002 Census Tract 310.02 Cobb County Georgia 33037 9787 3286 2853 MULTIPOLYGON (((-84.58191 3... ## 10 13089020200 Census Tract 202 DeKalb County Georgia 87889 2045 1704 187 MULTIPOLYGON (((-84.34896 3... ``` ]] --- ## What's inside geometry column .footnotesize[ .scroll-box-10[ ```r census %>% st_geometry() ``` ``` ## Geometry set for 519 features ## Geometry type: MULTIPOLYGON ## Dimension: XY ## Bounding box: xmin: -84.85071 ymin: 33.35246 xmax: -84.02371 ymax: 34.18629 ## Geodetic CRS: WGS 84 ## First 5 geometries: ``` ``` ## MULTIPOLYGON (((-84.37445 33.60568, -84.37159 3... ``` ``` ## MULTIPOLYGON (((-84.44869 33.57185, -84.43072 3... ``` ``` ## MULTIPOLYGON (((-84.42895 33.59714, -84.42687 3... ``` ``` ## MULTIPOLYGON (((-84.27173 33.56137, -84.26847 3... ``` ``` ## MULTIPOLYGON (((-84.65095 34.0392, -84.64618 34... ``` ] ] .scroll-box-10[.footnotesize[ ```r census %>% st_coordinates() ``` ``` ## X Y L1 L2 L3 ## [1,] -84.37445 33.60568 1 1 1 ## [2,] -84.37159 33.60570 1 1 1 ## [3,] -84.36760 33.60381 1 1 1 ## [4,] -84.36577 33.60348 1 1 1 ## [5,] -84.36407 33.60370 1 1 1 ## [6,] -84.36292 33.60384 1 1 1 ## [7,] -84.36180 33.60364 1 1 1 ## [8,] -84.35833 33.60305 1 1 1 ## [9,] -84.35526 33.60328 1 1 1 ## [10,] -84.35525 33.60044 1 1 1 ## [11,] -84.35434 33.59889 1 1 1 ## [12,] -84.35421 33.59881 1 1 1 ## [13,] -84.35341 33.59868 1 1 1 ## [14,] -84.35308 33.59858 1 1 1 ## [15,] -84.35175 33.59784 1 1 1 ## [16,] -84.35098 33.59671 1 1 1 ## [17,] -84.35189 33.59257 1 1 1 ## [18,] -84.35141 33.59153 1 1 1 ## [19,] -84.35110 33.59107 1 1 1 ## [20,] -84.35106 33.58568 1 1 1 ## [21,] -84.35253 33.58715 1 1 1 ## [22,] -84.35303 33.58750 1 1 1 ## [23,] -84.35438 33.58803 1 1 1 ## [24,] -84.35577 33.58829 1 1 1 ## [25,] -84.35834 33.58847 1 1 1 ## [26,] -84.35914 33.58863 1 1 1 ## [27,] -84.36085 33.58917 1 1 1 ## [28,] -84.36346 33.59025 1 1 1 ## [29,] -84.36369 33.59034 1 1 1 ## [30,] -84.36477 33.59064 1 1 1 ## [31,] -84.36590 33.59071 1 1 1 ## [32,] -84.36791 33.59079 1 1 1 ## [33,] -84.37007 33.59086 1 1 1 ## [34,] -84.37082 33.59086 1 1 1 ## [35,] -84.37184 33.59085 1 1 1 ## [36,] -84.37209 33.59306 1 1 1 ## [37,] -84.37328 33.59450 1 1 1 ## [38,] -84.37285 33.59692 1 1 1 ## [39,] -84.37411 33.60111 1 1 1 ## [40,] -84.37473 33.60367 1 1 1 ## [41,] -84.37445 33.60568 1 1 1 ## [42,] -84.44869 33.57185 1 1 2 ## [43,] -84.43072 33.57923 1 1 2 ## [44,] -84.42833 33.58026 1 1 2 ## [45,] -84.42755 33.57917 1 1 2 ## [46,] -84.42709 33.57849 1 1 2 ## [47,] -84.42632 33.57751 1 1 2 ## [48,] -84.42456 33.57752 1 1 2 ## [49,] -84.42128 33.57881 1 1 2 ## [50,] -84.41551 33.57883 1 1 2 ## [51,] -84.41203 33.57876 1 1 2 ## [52,] -84.41259 33.57265 1 1 2 ## [53,] -84.41316 33.56650 1 1 2 ## [54,] -84.41437 33.56505 1 1 2 ## [55,] -84.41686 33.56150 1 1 2 ## [56,] -84.42142 33.56197 1 1 2 ## [57,] -84.42525 33.56243 1 1 2 ## [58,] -84.42761 33.56375 1 1 2 ## [59,] -84.42930 33.56427 1 1 2 ## [60,] -84.43312 33.56615 1 1 2 ## [61,] -84.43481 33.56669 1 1 2 ## [62,] -84.44808 33.56701 1 1 2 ## [63,] -84.44783 33.56913 1 1 2 ## [64,] -84.44869 33.57185 1 1 2 ## [65,] -84.42895 33.59714 1 1 3 ## [66,] -84.42687 33.59714 1 1 3 ## [67,] -84.42340 33.59613 1 1 3 ## [68,] -84.41631 33.59578 1 1 3 ## [69,] -84.41435 33.59494 1 1 3 ## [70,] -84.41241 33.59462 1 1 3 ## [71,] -84.40846 33.59408 1 1 3 ## [72,] -84.40795 33.59392 1 1 3 ## [73,] -84.40958 33.58856 1 1 3 ## [74,] -84.41003 33.58716 1 1 3 ## [75,] -84.41686 33.58728 1 1 3 ## [76,] -84.42172 33.58696 1 1 3 ## [77,] -84.42330 33.58817 1 1 3 ## [78,] -84.42637 33.58989 1 1 3 ## [79,] -84.42897 33.59293 1 1 3 ## [80,] -84.42895 33.59714 1 1 3 ## [81,] -84.27173 33.56137 1 1 4 ## [82,] -84.26847 33.56358 1 1 4 ## [83,] -84.26739 33.56601 1 1 4 ## [84,] -84.26438 33.56779 1 1 4 ## [85,] -84.26339 33.56933 1 1 4 ## [86,] -84.26347 33.56050 1 1 4 ## [87,] -84.26385 33.55521 1 1 4 ## [88,] -84.26392 33.55382 1 1 4 ## [89,] -84.26370 33.54892 1 1 4 ## [90,] -84.26864 33.54896 1 1 4 ## [91,] -84.26922 33.54896 1 1 4 ## [92,] -84.26894 33.55315 1 1 4 ## [93,] -84.26942 33.55550 1 1 4 ## [94,] -84.27173 33.56137 1 1 4 ## [95,] -84.65095 34.03920 1 1 5 ## [96,] -84.64618 34.03719 1 1 5 ## [97,] -84.64582 34.04214 1 1 5 ## [98,] -84.64575 34.04221 1 1 5 ## [99,] -84.64511 34.04610 1 1 5 ## [100,] -84.64338 34.04656 1 1 5 ## [101,] -84.64382 34.04688 1 1 5 ## [102,] -84.64428 34.04758 1 1 5 ## [103,] -84.64546 34.05095 1 1 5 ## [104,] -84.64267 34.05224 1 1 5 ## [105,] -84.64213 34.05198 1 1 5 ## [106,] -84.64142 34.05160 1 1 5 ## [107,] -84.63989 34.05128 1 1 5 ## [108,] -84.63795 34.05038 1 1 5 ## [109,] -84.63130 34.04373 1 1 5 ## [110,] -84.62708 34.04133 1 1 5 ## [111,] -84.62750 34.03338 1 1 5 ## [112,] -84.62791 34.03179 1 1 5 ## [113,] -84.62689 34.02943 1 1 5 ## [114,] -84.62513 34.02613 1 1 5 ## [115,] -84.62521 34.02209 1 1 5 ## [116,] -84.63560 34.02634 1 1 5 ## [117,] -84.63753 34.02694 1 1 5 ## [118,] -84.63850 34.02715 1 1 5 ## [119,] -84.65000 34.02942 1 1 5 ## [120,] -84.65000 34.02953 1 1 5 ## [121,] -84.64999 34.03154 1 1 5 ## [122,] -84.65100 34.03390 1 1 5 ## [123,] -84.65095 34.03920 1 1 5 ## [124,] -84.70304 33.94043 1 1 6 ## [125,] -84.69364 33.94382 1 1 6 ## [126,] -84.68672 33.94481 1 1 6 ## [127,] -84.68022 33.94530 1 1 6 ## [128,] -84.67741 33.94588 1 1 6 ## [129,] -84.67334 33.94739 1 1 6 ## [130,] -84.66163 33.95244 1 1 6 ## [131,] -84.66240 33.94951 1 1 6 ## [132,] -84.66392 33.94829 1 1 6 ## [133,] -84.66407 33.94655 1 1 6 ## [134,] -84.66369 33.94177 1 1 6 ## [135,] -84.66396 33.93898 1 1 6 ## [136,] -84.66180 33.93480 1 1 6 ## [137,] -84.66200 33.92422 1 1 6 ## [138,] -84.66204 33.92402 1 1 6 ## [139,] -84.66385 33.92319 1 1 6 ## [140,] -84.66544 33.92132 1 1 6 ## [141,] -84.66741 33.91434 1 1 6 ## [142,] -84.66654 33.91009 1 1 6 ## [143,] -84.66660 33.90654 1 1 6 ## [144,] -84.66703 33.90223 1 1 6 ## [145,] -84.67043 33.90235 1 1 6 ## [146,] -84.67264 33.90251 1 1 6 ## [147,] -84.67335 33.90483 1 1 6 ## [148,] -84.67449 33.90565 1 1 6 ## [149,] -84.67945 33.90777 1 1 6 ## [150,] -84.68273 33.90952 1 1 6 ## [151,] -84.68522 33.91131 1 1 6 ## [152,] -84.69111 33.91695 1 1 6 ## [153,] -84.69224 33.91861 1 1 6 ## [154,] -84.69939 33.93236 1 1 6 ## [155,] -84.70185 33.93908 1 1 6 ## [156,] -84.70304 33.94043 1 1 6 ## [157,] -84.55734 34.05793 1 1 7 ## [158,] -84.54304 34.05750 1 1 7 ## [159,] -84.53826 34.05741 1 1 7 ## [160,] -84.53634 34.05657 1 1 7 ## [161,] -84.52867 34.04992 1 1 7 ## [162,] -84.52733 34.04928 1 1 7 ## [163,] -84.52952 34.04580 1 1 7 ## [164,] -84.52957 34.04401 1 1 7 ## [165,] -84.52827 34.03844 1 1 7 ## [166,] -84.52766 34.03547 1 1 7 ## [167,] -84.52777 34.03138 1 1 7 ## [168,] -84.52786 34.02977 1 1 7 ## [169,] -84.52849 34.02564 1 1 7 ## [170,] -84.52845 34.02068 1 1 7 ## [171,] -84.52833 34.01996 1 1 7 ## [172,] -84.53142 34.01908 1 1 7 ## [173,] -84.53454 34.01893 1 1 7 ## [174,] -84.53737 34.01578 1 1 7 ## [175,] -84.53988 34.01420 1 1 7 ## [176,] -84.54317 34.01378 1 1 7 ## [177,] -84.54841 34.01382 1 1 7 ## [178,] -84.54830 34.01967 1 1 7 ## [179,] -84.54795 34.02115 1 1 7 ## [180,] -84.54780 34.02199 1 1 7 ## [181,] -84.54815 34.02339 1 1 7 ## [182,] -84.54949 34.02608 1 1 7 ## [183,] -84.55022 34.03402 1 1 7 ## [184,] -84.55078 34.03990 1 1 7 ## [185,] -84.55238 34.04310 1 1 7 ## [186,] -84.55557 34.04655 1 1 7 ## [187,] -84.55608 34.04761 1 1 7 ## [188,] -84.55618 34.05213 1 1 7 ## [189,] -84.55653 34.05543 1 1 7 ## [190,] -84.55773 34.05790 1 1 7 ## [191,] -84.55734 34.05793 1 1 7 ## [192,] -84.51010 33.93943 1 1 8 ## [193,] -84.50777 33.93943 1 1 8 ## [194,] -84.50712 33.93797 1 1 8 ## [195,] -84.50447 33.93605 1 1 8 ## [196,] -84.50277 33.93392 1 1 8 ## [197,] -84.50122 33.93313 1 1 8 ## [198,] -84.50051 33.93280 1 1 8 ## [199,] -84.50042 33.93277 1 1 8 ## [200,] -84.49810 33.93185 1 1 8 ## [ reached getOption("max.print") -- omitted 17081 rows ] ``` ]] --- ## Interactive Mapping - tmap package ```r library(tmap) tmap_mode('view') ``` ``` ## tmap mode set to interactive viewing ``` ```r tm_shape(census) + tm_polygons(col="hhincome") ```
--- ## Data manipulation .font90[ #### sf objects ARE data frames. .font60[.gray[(Note that sfc, sfg, etc. are not dataframes.)]] Most verbs from dplyr packages work with sf objects. .scroll-box-18[ ```r census_group <- census %>% select(county, hhincome) %>% filter(county %in% c("Clayton County", "Fulton County")) %>% group_by(county) %>% summarise(count = n(), med.income = median(hhincome, na.rm = TRUE)) census_group ``` ``` ## Simple feature collection with 2 features and 3 fields ## Geometry type: POLYGON ## Dimension: XY ## Bounding box: xmin: -84.85071 ymin: 33.35246 xmax: -84.09769 ymax: 34.18629 ## Geodetic CRS: WGS 84 ## # A tibble: 2 × 4 ## county count med.income geometry ## <chr> <int> <dbl> <POLYGON [°]> ## 1 Clayton County 50 43321 ((-84.45893 33.56541, -84.45892 33.55976, -84.45858 33.55584, -84.45866 33.55093, -84.448... ## 2 Fulton County 204 56944 ((-84.17983 34.008, -84.17906 34.0091, -84.17832 34.0123, -84.17664 34.01627, -84.17529 3... ``` ]] --- The **group_by() + summarise()** combination modifies the geometry too. ```r tm_shape(census_group) + tm_polygons(col = "med.income") ```
--- ## Subsetting by location When you have two or more sf objects, you can select a subset. Let's do a demo with Census Tracts and a 10-KM buffer from the Five Points Station. .pull-left[ .footnotesize[ ```r # Visualization leaflet() %>% addProviderTiles(providers$CartoDB.Positron) %>% * addPolygons(data = census, weight=1, color="grey") %>% * addPolygons(data = five_point, weight=1) %>% addLegend("bottomright", colors = c("blue", "grey"), labels = c("five_point", "poly"), title="Legend") ``` It's okay if you don't understand this code. We will learn about **leaflet** package later. ]] .pull-right[ .script-size[
]] --- ## Subsetting The magic syntax is: .center[.red[sf-obj-1][.blue[sf-obj-2], , .green[operation]]] .small[ ```r # Subsetting census.subset <- census[five_point, , op = st_intersects] ``` ] .scriptsize[ ```r # Visualization tm_shape(census.subset) + tm_polygons() + tm_shape(five_point) + tm_borders(col = "blue") ```
] --- ## Predicates .font80[ In GIS, "**predicates are Boolean functions** that return TRUE if a comparison meets the functions criteria; otherwise, they return FALSE (f), to determine if a specific relationship exists between a pair of geometries." ([ESRI](https://help.arcgis.com/en/geodatabase/10.0/sdk/arcsde/concepts/geometry/shapes/spatial_relations/predicates.htm)). Try **?st_intersects** in your console. ] <img src="predicates.jpg" width="80%" /> --- .footnotesize[ ```r a <- leaflet() %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons(data = pred_is_within_distance, weight=1,color="grey") %>% addPolygons(data = five_point, color="black") %>% addControl(html="<b>st_is_with_distance (10k)</b>", "topright") b <- leaflet() %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolygons(data = pred_disjoint, weight=1,color="grey") %>% addPolygons(data = five_point, color="black") %>% addControl(html="<b>st_disjoint</b>", "topright") leafsync::sync(a,b) ```
] --- ## Predicates <!-- --> --- ## Buffer, intersection, etc. Usual GIS operations work the same way in R as they do in ArcMap, for example. .footnotesize[ ```r # Buffer (5K) + intersection buf <- marta %>% st_buffer(dist = units::set_units(5000, "m")) buf_int_census <- buf %>% st_intersection(census) ``` ``` ## Warning: attribute variables are assumed to be spatially constant throughout all geometries ``` ```r tm_shape(census) + tm_polygons(col = "gray") + tm_shape(buf_int_census) + tm_polygons(col = "red") + tm_shape(marta) + tm_lines(col = "blue", lwd = 3) ```
] --- class: inverse, center, middle #If we have time left.. --- ## Visualizing linestring and points Let's visualize the first LINESTRING (i.e., Blue Line of MARTA) and the points that makes the LINESTRING. .pull-left[ #### Extracting points .scroll-box-20[ .footnotesize[ ```r marta.point <- marta %>% # Extract the first row slice(2) %>% # Turn it into coordinates st_coordinates() %>% # Turn it into df as.data.frame() %>% # df to sf st_as_sf(coords = c("X", "Y"), crs = 4326) marta.point ``` ``` ## Simple feature collection with 52 features and 2 fields ## Geometry type: POINT ## Dimension: XY ## Bounding box: xmin: -84.44894 ymin: 33.63922 xmax: -84.28036 ymax: 33.90298 ## Geodetic CRS: WGS 84 ## First 10 features: ## L1 L2 geometry ## 1 1 1 POINT (-84.44628 33.63922) ## 2 1 1 POINT (-84.44894 33.65525) ## 3 2 1 POINT (-84.44894 33.65525) ## 4 2 1 POINT (-84.44017 33.67876) ## 5 3 1 POINT (-84.44017 33.67876) ## 6 3 1 POINT (-84.44008 33.6852) ## 7 4 1 POINT (-84.44008 33.6852) ## 8 4 1 POINT (-84.42975 33.69562) ## 9 5 1 POINT (-84.42975 33.69562) ## 10 5 1 POINT (-84.42379 33.72086) ``` ] ] ] .pull-right[ #### Extracing line .scroll-box-20[ .footnotesize[ ```r marta.line <- marta %>% # Extract the first row slice(2) marta.line ``` ``` ## Simple feature collection with 1 feature and 2 fields ## Geometry type: MULTILINESTRING ## Dimension: XY ## Bounding box: xmin: -84.44894 ymin: 33.63922 xmax: -84.28036 ymax: 33.90298 ## Geodetic CRS: WGS 84 ## route_long OBJECTID geometry ## 1 GOLD-Northeast Doraville Line 1552.5 MULTILINESTRING ((-84.44628... ``` ] ] ] --- ## Plotting them together ```r library(leaflet) leaflet() %>% addProviderTiles(providers$CartoDB.Positron) %>% # Background map addPolylines(data = marta %>% slice(c(1,3)), color="grey", weight=1) %>% # Adding non-Blue lines * addPolylines(data = marta.line) %>% # Adding marta.line * addCircles(data = marta.point) # Adding marta.point ```
--- ## Why this is useful Being able to parse the geometry at point level allows us to do any manipulation we want to do. For example, we want to cut the line into two at near Five Points Station (33.7539, -84.3916). .footnotesize[ .scroll-box-14[ ```r # Creating POINT feature for five point station five_point_station <- st_point(c(-84.3916, 33.7539)) %>% st_sfc(crs = 4326) # Calculate pairwise distance to all points that define Gold line dist.mat <- marta.point %>% st_distance(five_point_station) # Which row contains the closest point? closest.pnt <- which.min(dist.mat) # Parse marta.line and create a cut-off line marta.line.cut <- marta.line %>% st_coordinates() %>% as.data.frame() %>% mutate(break.point = c(rep("South",closest.pnt), rep("North",nrow(.)-closest.pnt))) %>% st_as_sf(coords=c("X", "Y"), crs = 4326) %>% group_by(break.point, L1) %>% summarise(do_union=FALSE) %>% st_cast("LINESTRING") %>% group_by(break.point) %>% summarise(do_union=FALSE) %>% st_cast("MULTILINESTRING") ``` ``` ## `summarise()` has grouped output by 'break.point'. You can override using the `.groups` argument. ``` ```r marta.line.cut ``` ``` ## Simple feature collection with 2 features and 1 field ## Geometry type: MULTILINESTRING ## Dimension: XY ## Bounding box: xmin: -84.44894 ymin: 33.63922 xmax: -84.28036 ymax: 33.90298 ## Geodetic CRS: WGS 84 ## # A tibble: 2 × 2 ## break.point geometry ## <chr> <MULTILINESTRING [°]> ## 1 North ((-84.38783 33.75747, -84.38608 33.77731), (-84.38608 33.77731, -84.38781 33.7917... ## 2 South ((-84.44628 33.63922, -84.44894 33.65525), (-84.44894 33.65525, -84.44017 33.6787... ``` ] ] .center[ ### No need to understand this code! I just wanted to demonstrate what can be done. ] --- .scriptsize[ ```r pal <- colorFactor(palette = c("red", "blue"), domain = marta.line.cut$break.point) leaflet() %>% addProviderTiles(providers$CartoDB.Positron) %>% addPolylines(data = marta.line.cut, color = ~pal(break.point)) %>% addPolylines(data = marta %>% slice(1,3), color="grey", weight=1) %>% addCircles(data=five_point_station, color="black", radius=300) %>% addLegend("bottomright", pal = pal, values = marta.line.cut$break.point) ```
]